DFT + U Study of Uranium Dioxide and Plutonium Dioxide with Occupation Matrix Control

DFT + U with occupation matrix control (OMC) is applied to study computationally bulk UO2 and PuO2, the latter for the first time. Using the PBESol functional in conjunction with OMC locates AFM and NM ground states for UO2 and PuO2, respectively, in agreement with experimental findings. By simulating the lattice parameter, magnetic moment, band gap, and densities of states, U = 4.0 eV is recommended for AFM UO2, yielding data close to experiments for all considered properties. U = 4.5 and 4.0 eV are recommended for NM and AFM PuO2, respectively, though much larger U values (c. 10 eV) are required to yield the most recently reported PuO2 band gap. For both oxides, several excited states have similar properties to the ground state, reinforcing the need to employ OMC wherever possible.


■ INTRODUCTION
In the actinide series of elements, many of the chemical and physical properties display a turning point at plutonium. This includes the change from the more covalent early actinides to the more ionic mid and later elements, and the extensive range of oxidation states exhibited by the early elements diminishes significantly after Pu. 1 The actinide dioxides, which are the subject of this work, change from Mott−Hubbard insulators to charge transfer insulators at PuO 2 . 2−4 UO 2 and NpO 2 have antiferromagnetic (AFM) ground states 5−7 while a nonmagnetic (NM) ground state is found for PuO 2 . 8,9 As PuO 2 is a product of the recycling of spent UO 2 nuclear fuel, detailed understanding of PuO 2 is clearly essential not just at a fundamental level but also to inform its safe current and longterm storage.
Due to the high radioactivity of PuO 2 , experiments are very challenging, and hence theoretical simulations play a particularly valuable role in its study. Density functional theory (DFT) with a Hubbard U correction is widely used 4,10,11 as it gives reasonable predictions at the lowest computational costs. However, the ability of the DFT + U approach to correctly identify the NM magnetic ground state of PuO 2 remains an issue. The NM ground state has been established by various experiments over a wide temperature range (4−1000 K), including by inelastic neutron scattering and nuclear magnetic resonance. 8,9,12−14 By contrast, previous DFT + U simulations have predicted an AFM ground state for PuO 2 ; 11,15,16 although this does not match with experiments, many of the other calculated properties of AFM PuO 2 do agree quite well. There are also some other theoretical works that adopt the experimentally indicated NM state. 4,17,18 The inconsistency between experiments and DFT + U simulation over the correct magnetic ground state of PuO 2 requires further study.
The majority of actinide compounds are open-shell and frequently feature several unpaired electrons in the seven valence 5f orbitals. There are typically many different ways in which the actinide f orbitals may be populated, and use of the Hubbard U parameter in DFT calculations can lead to the location of excited states arising from those electronic configurations. 19,20 It is unclear whether the previous computational reports of an AFM ground state for PuO 2 arise from a fundamental inability of DFT + U to locate the correct ground state or if they have become trapped in higher energy states. To address this question, we here consider all possible filling patterns of the 4 electrons in the 7 5f orbitals of Pu(IV), using the occupation matrix control (OMC) approach. 21 To the best of our knowledge, there is no such work so far. Wang and Konashi considered all possible occupation matrices (OMs) for AFM Pu (5f 5 ) in PuO 2 , 22 though not Pu (5f 4 ) in which we are interested, and found that there are low-lying excited states. Dorado et al. considered all possible OMs for the 2 electrons in the 7 5f orbitals of AFM UO 2 , 23 finding that the highest energy OM state is about 3.5 eV above the lowest one.
In this work, we first consider all possible OMs for AFM, FM, and NM UO 2 . Although previous DFT + U simulations agree with experiments over the AFM ground state for UO 2 , 24 it is worth considering all possible OMs for FM and NM as well as AFM to be sure of the correct computed ground state for UO 2 , and UO 2 provides a good test of the DFT + U with the OMC method as there are more data available on UO 2 . We then consider all possible OMs for AFM, FM, and NM PuO 2 , in order to establish the correct theoretical ground state. Through exploration of the effect of the choice of Hubbard U on a range of computed properties (lattice parameter, band gap, magnetic moment, and density of states), we also aim to provide recommendations as to the best values of U to employ in DFT + U + OMC studies of UO 2 and PuO 2 .

■ COMPUTATIONAL DETAILS
All calculations were performed using density functional theory (DFT), as implemented in the Vienna ab initio simulation package (VASP), version 5.4.1. 25−28 The generalized gradient approximation functional of Perdew, Burke, and Ernzerhof, revised for solids (PBESol), was used, 29 with a Hubbard U correction for the 5f electrons. 30 A wide range of U values (0.0−7.0 eV) was considered to establish the most suitable values for the simulation of UO 2 and PuO 2 bulk (Figure 1a). Plane wave basis sets and projector augmented wave pseudopotentials were used to describe the ions. 31 Plane wave cutoff energy and k mesh sizes were tested for UO 2 bulk with the lattice parameter fixed at 5.470 Å (the experimental values are 5.470−5.473 Å); 32−34 Figure 1b shows that a plane wave cutoff of 500 eV and gamma-centered 5 × 5 × 5 Monkhorst−Pack grid for the Brillouin zone are sufficient. 35 Therefore, a 650 eV (1.3 × 500 eV) cutoff energy (to weaken the influence of Pulay stress) and gamma-centered 5 × 5 × 5 Monkhorst−Pack grid were used for all calculations in this work. The iteration threshold for electronic and ionic convergence was set to 1 × 10 −5 and 1 × 10 −2 eV, respectively. Antiferromagnetic (AFM), ferromagnetic (FM), and nonmagnetic (NM) states were considered with 1 k colinear ordering (along the c direction, Figure 1a) for both UO 2 and PuO 2 . 1 k colinear ordering was chosen over 3 k non-colinear ordering because, although AnO 2 exhibit non-collinear magnetic behavior 36 in which the magnetic moments of the ions have contributions in more than one direction, 1 k ordering is much more computationally tractable than 2 k or 3 k ordering. Furthermore, most previous computational work also uses 1 k ordering, so our using it facilitates more direct comparison, and we here in part aim to find a theoretical approach that gives accurate simulation of UO 2 and PuO 2 at manageable computational cost. As the 1 k colinear ordering is used, we have only Type G magnetic arrangements (along the c direction). U/Pu labeled 1 and 2 in Figure 1a are set to spin up and 3 and 4 are set to spin down for all AFM calculations in this work.
Occupation matrix control (OMC), developed by Dorado et al. 23,37 and incorporated into VASP by Allen and Watson, 21 was used to explore all possible OMs. Only the diagonal elements were set to non-zero values for initial OMs: i k j j j j j j j j j j j j j j j j j j j j j j j j j j j j j y { z z z z z z z z z z z z z z z z z z z z z z z z z z z z z where n is either 0.0 or 1.0. U and Pu in their dioxide bulk have 2 and 4 5f electrons, respectively; these are unpaired in the FM and AFM states and paired in the NM state. Therefore, we studied C 7 2 = 21 OMs for FM and AFM UO 2 , C 7 1 = 7 OMs for NM UO 2 , C 7 4 = 35 OMs for FM and AFM PuO 2 , and C 7 2 = 21 OMs for NM PuO 2 ; all OMs are listed in the Supporting Information (Tables S1 and S2). Although some of the electronic configurations defined by the OMs are degenerate, we decided to follow the approach of Dorado et al., who used DFT + U + OMC to investigate bulk UO 2 23 and study all the OMs as the imposition of U can decrease the degeneracy of the f orbitals. The initially imposed OMs remain unchanged during the self-consistent field calculations.
All of the data presented in the main text were obtained with the PBESol functional. We also performed all the calculations using the PBE functional and found that PBESol predicts better lattice parameters. The PBE data are collected in Figures S5−S9 and Table S3 in the Supporting Information.
■ RESULTS AND DISCUSSION Uranium Dioxide (UO 2 ). We now explore the energies of the AFM, FM, and NM states of UO 2 with the OMs as described in the Computational Details. As the Hubbard U has an influence on the localization and energy of the f orbitals, a wide range of U values is considered. The pure PBEsol method (U = 0.0 eV) predicts the same energy for all solutions of AFM, FM and NM UO 2 ( Figure 2). Furthermore, checking the magnetic moment on each U atom indicates that all AFM and NM states optimize to FM states. Introduction of a non- The Journal of Physical Chemistry C pubs.acs.org/JPCC Article zero U value breaks the degeneracy of the f orbitals obtained from the pure DFT simulation, the data for U = 1.0 to 7.0 eV in Figure 2 are no longer linear, and stable AFM and NM solutions are found. The effect of the Hubbard U becomes increasingly clear with increasing U value. Energy differences between different solutions are amplified at high U, as wider energy ranges of solutions are found. Small U values (<2 eV) are not enough to overcome the drawback of pure DFT as some AFM and NM states still optimize to FM states, but larger U values (≥ 2 eV) predict the same OM for the most stable solution. For AFM and FM, the OM has two unpaired electrons occupying the f 1 and f 3 orbitals (the 20th OM listed in Table S1), while for the NM state, the OM has pair of electrons in the f 1 orbital (the 5th OM listed in Table S1). In the following discussion, we focus on the most stable AFM, FM, and NM OM solutions, and refer to them as the AFM, FM, and NM states for simplicity.
As noted above, pure DFT and low U values (<2.0 eV) predict an FM ground state, regardless of the initially chosen magnetic state; the energies obtained with U ≥ 2.0 eV are compared for AFM, FM, and NM states in Figure 3. As shown in Figure 3b, an NM ground state is found for UO 2 with 2.0 ≤ U < 3.0 eV; when U ≥ 3.0 eV, the AFM ground state appears. Previous experiments have confirmed that the ground state of UO 2 is AFM, so the chosen U value should be not smaller than 3.0 eV to give reasonable prediction. As the ground state is predicted differently with low and high U values, we suggest caution when employing the U-ramping method, in which U is scanned from zero to the desired value in small steps, while reading the previous step's result.
We did not find any distortion of cubic UO 2 bulk, i.e., 1 k AFM UO 2 keeps the Fm-3m crystal symmetry. The lattice parameter gradually increases with the U value, and 3 eV ≤ U ≤ 4 eV yields lattice parameters close to the experimental values ( Figure 4a). The magnetic moment on the U ions also increases with increasing U (Figure 4b); imposing U on the f orbitals leads to more localized f electrons. AFM and FM UO 2 simulations slightly overestimate the magnetic moment on U ions, but the differences are small, and AFM UO 2 is predicted to have a magnetic moment closer to the experimental value. 38 The band gap of UO 2 increases with U in a similar way for the AFM, FM, and NM states. Across the range of U studied, NM UO 2 has the largest band gap and FM UO 2 has the smallest band gap, with the band gap of AFM UO 2 being slightly lower than for NM. Figure 4c suggests that 4.0 eV ≤ U ≤ 5.0 eV gives good agreement with experiments for the band gap of AFM UO 2 . 39,40 As well as the band gap, the density of states (DOS) is an important electronic property against which to evaluate the simulation. As AFM is the experimentally and theoretically reported ground state, only the DOS of AFM UO 2 is discussed here. We plot the DOS for AFM UO 2 with U values ranging from 2.0 to 7.0 eV in steps of 1.0 eV ( Figure 5). We are mainly interested in three bands in the DOS: the valence band, the conduction band, and the second band below the Fermi level. Increasing the Hubbard U localizes and stabilizes the f electrons. The valence band, which is mainly U 5f, moves downward in the DOS. This increases the band gap, and the gap between the valence band and the second band (which is mainly of O 2p character) reduces, such that the two bands become mixed at high U (6.0 and 7.0 eV). As the Fermi level is fixed at 0.0 eV, the conduction band and the second band under the Fermi level move upward with increasing U value.
X-ray absorption data indicate that the valence band is mainly U 5f, 41 so we can exclude U ≥ 6.0 eV. Previous experiments also show that the second peak under the Fermi level consists of O 2p states at around −4.0 eV, 41,42 so we can also exclude U = 5.0 eV. Though U = 2.0 eV predicts a reasonable position for the second peak under the Fermi level, the band gap is much smaller than the experimental value and it predicts an NM ground state. Therefore, 2.0 eV < U < 5.0 eV is the preferred range to obtain reasonable DOS for UO 2 bulk simulation.
Overall, we find that 3.0 eV ≤ U ≤ 5.0 eV gives the best balance of agreement with experimental data over a range of properties, with the best value being 4.0 eV. Some properties of AFM UO 2 calculated with the PBESol + U (U = 4.0 eV) + OMC method are summarized and compared with experimental values in Table 1.     Table S2), while for the NM states, the most stable solution has two pairs of electrons in the f 1 and f 3 orbitals (the 20th OM in Table S2). In the following discussion, we focus on the most stable AFM, FM, and NM solutions and refer to them as the AFM, FM, and NM states for simplicity. Properties of the AFM, FM, and NM states are compared in Figure 6. When U ≥ 2.0 eV, an NM ground state is found, in agreement with previous experiments in the temperature range 4−1000 K. 8,14 However, the energy difference between AFM and NM is small (Figure 6b), increasing from U = 2.0 to 4.0 eV and then decreasing again, with the largest difference of −0.18 eV. Previous DFT + U simulations have found an AFM ground state for PuO 2 , although given the small energy differences with NM states, it may be that DFT + U simulation of PuO 2 without OMC can become trapped in an AFM state.
Similar to UO 2 , PuO 2 bulk remains face-centered cubic after optimization and the optimized lattice parameter increases with the U value (Figure 6c   The Journal of Physical Chemistry C pubs.acs.org/JPCC Article the same U value. This is by contrast to UO 2 , where although the largest lattice parameter was also obtained for the FM state, it is close to the AFM and NM lattice parameters for a given U value (Figure 4a). This may well be reminiscent of previous work on the paramagnetic to ferromagnetic transition of La(Fe x Si 1−x ) 13 , 47 where a larger volume change is observed for x = 0.88 (where Fe has a larger magnetic moment) than for x = 0.86 (where Fe has a smaller magnetic moment); as Pu 4+ has a larger magnetic moment than U 4+ , the difference between the lattice parameter of FM and AFM/NM for PuO 2 is larger than for UO 2 . PuO 2 is found to be nonmagnetic by experiment, and hence the magnetic moment of Pu should be zero, and experiments indeed find only a very small nuclear magnetic moment for Pu 4+ (about 0.15 μ N ), 14 which may arise from coupling between the singlet Γ 1 ground state and (an) excited state(s). The excited state could be solely the triplet Γ 4 state, at an energy of about 0.120 eV, or two or more states in the energy region 0.110−0.140 eV, but spectral resolution is insufficient to be certain. 9 We here find that AFM PuO 2 is higher in energy than NM PuO 2 by less than 0.18 eV, so ground-state PuO 2 could have a small contribution from the AFM state; highresolution neutron spectroscopy would be helpful here. The Pu magnetic moments of AFM and FM PuO 2 against the U value are given in Figure S2; due to more localized f states with increasing U value, the magnetic moment of Pu increases with U.
Pure PBESol and PBESol + U with a small U value (< 3.0 eV) predict PuO 2 to be metallic; with increasing U, PuO 2 becomes a semiconductor with a band gap that increases with U. As with UO 2 , NM PuO 2 has the largest band gap, and FM has the smallest. The valence band of PuO 2 is o mixed f and p character, with a higher contribution of the former, 50−52 so U ≥ 5.0 eV values are not good choices as they predict a higher (or equal) contribution of p states than f states for the valence band. Previous X-ray photoelectron spectra have shown that the valence band is split into a mainly Pu 5f-contributed state (at higher energy) and a mainly O 2p-contributed state (at lower energy), 50 which is also supported by previous theoretical simulations. Previous UPS studies showed that the contribution of Pu 5f states is centered at around 2 eV below the Fermi level, in agreement with XPS data. 51 XPS also shows that the mainly O 2p band, which covers a wide range below the Fermi level, extends to −10 eV, 50 while it ends at around −8 eV from UPS data. 51 The center of the Pu 5f states moves downward with increasing U; for NM PuO 2 simulated with U = 2.0 eV, the 5f states are centered at around −1.5 eV, at around −2.0 eV with U = 4.0 eV, matching the experiment well, and at around −2.5 eV with U = 5.0 eV, which is also reasonable. The O 2p states end at around −7 eV for NM PuO 2 with U = 2.0 eV, which is the closest value to experimental data (about −8 eV), as the second band moves upward with increasing U, ending at around −6 eV for NM PuO 2 with U = 5.0 eV. In general, 4.0 eV ≤ U < 5.0 eV gives reasonable prediction for the position of the Pu 5f states in the valence band (Figure 7). We have also The Journal of Physical Chemistry C pubs.acs.org/JPCC Article studied NM PuO 2 bulk with U = 4.5 eV; the DOS is given in Figure 8, which meets the experimentally reported features well, such as the composition and center of the valence band. Previous works also suggest that there is a peak of the O p state character on the left shoulder of the valence band, a feature which can be seen in our DOS for NM PuO 2 simulated with U = 4.0 and 4.5 eV. Overall, 4.0 eV ≤ U ≤ 4.5 eV is good for NM PuO 2 bulk simulation with U = 4.5 eV being the best, as this value predicts a larger band gap than U = 4.0 eV. We also studied the DOS of AFM PuO 2 as it is the most studied state of PuO 2 and may have a contribution to the ground state. The DOS of AFM PuO 2 ( Figure 9) is similar to those of NM PuO 2 , with some minor differences. The DOS of AFM PuO 2 shown here is also similar to previous theoretical simulations. Therefore, although most previous DFT + U works without OMC study the AFM state of PuO 2 , they still obtain results similar to experiments, i.e., AFM PuO 2 is a good approximation to NM PuO 2 in the simulation of certain properties. To obtain good DOS for AFM PuO 2 , the U value should be smaller than 5 eV as U ≥ 5.0 eV predicts an O 2p state-dominated valence band. The DOS of PuO 2 calculated with U = 4.5 eV is given in Figure S3 and predicts almost the same contribution of p and f states to the valence band, so U ≤ 4.5 eV is needed to predict an f state-dominated valence band, while U < 3.0 eV predicts AFM PuO 2 as metallic. Overall, for AFM PuO 2 , 4.0 eV ≤ U ≤ 4.5 eV with U = 4.0 eV is the best as it predicts more reasonable DOS.
In summary, first, to obtain an NM ground state for PuO 2 , the chosen U value must be larger than 2 eV; second, to obtain reasonable lattice parameters, 4.0 eV ≤ U ≤ 5.0 eV is suggested; third, 5.0 eV ≤ U ≤ 10.0 eV is needed to reproduce the range of experimentally reported band gaps; fourth, to give reasonable DOS, 4.0 eV ≤ U ≤ 5.0 eV is suggested. Hence, no single U value can simultaneously reproduce the band gap and the other three pieces of experimental data. U = 10.0 predicts a band gap of 4.10 eV, which is the latest experimental datum, but it also predicts totally wrong DOS ( Figure S4). Overall, therefore, we suggest 4.0 eV ≤ U ≤ 5.0 eV with 4.5 and 4.0 eV being the best values for the simulation of NM and AFM PuO 2 bulk, respectively. Some properties of NM and AFM PuO 2 calculated with PBESol + U (4.5/4.0 eV) + OMC are listed in Table 1.
Dependence of Energy, Lattice Parameter, and Band Gap on the Occupation Matrix. With the ground state OMs and the effect of the Hubbard U on the properties of those ground states in hand for both UO 2 and PuO 2 , we here provide an insight into how the solutions calculated from different initial OMs affect the energy, lattice parameter, and band gap. Table 2 presents these data for the solutions arising from the different initial OMs of AFM UO 2 and NM PuO 2 , using U values of 4.0 and 4.5 eV, respectively.
We find the same OM for ground-state AFM UO 2 as does previous work, 23 and the energy range spanned by the solutions is about 9 eV. This is significantly larger than in the previous study, where the highest energy solution is only 3.45 eV higher than the ground state, although we note that previous workers report some non-converged states, which have the same OMs as the high-energy solutions (Table 2, in red) we find here. These high-energy solutions have a slightly larger lattice parameter (∼5.53 Å) than the ground state (5.479 Å) and much smaller band gaps (∼1.0 eV vs 1.99 eV). We expect that, due to their high energy, optimizations without initial OMC have a good chance of avoiding these solutions.
There are 6 medium-energy solutions of AFM UO 2 with relative energies of c. 1 to 3 eV (Table 2, in orange) and a further 8 low-energy solutions with relative energies of <1 eV    Table 2, in orange), and 4 (Table 2, in blue) high-, medium-, and lowenergy solutions, respectively. The high-energy solutions have relative energies of c. 7−28 eV, i.e., they are very unstable, while the medium-and low-energy solutions have relative energies of c. 1−4 eV and < 1 eV, respectively. Most of the medium-and low-energy solutions of NM PuO 2 have significantly larger or smaller band gaps than the ground state; there is only one medium-energy solution ([0002200]) and one low-energy solution ([2020000]) with similar lattice parameters and band gaps to the ground state.
Overall, for both oxides, there are some solutions with similar energies and lattice parameters to the ground state but with very different band gaps. However, other states have similar energies, lattice parameters, and band gaps to the ground state, and hence it is very hard to distinguish them from the ground state without OMC.

■ CONCLUSIONS
In this work, we have studied bulk UO 2 and PuO 2 with PBESol + U (0−7 eV) + OMC. By calculating the energies of all possible solutions with different initially imposed OMs of 1 k AFM, FM and NM UO 2 and PuO 2 , PBESol + U + OMC simulation predicts AFM and NM ground states for UO 2 and PuO 2 , respectively. Our UO 2 ground state is in agreement with previous experimental results and theoretical simulations. For PuO 2 , we show for the first time that PBESol + U + OMC correctly reproduces the experimentally reported NM ground state. We have also considered a wide range of U in order to find the best value for theoretical simulation. The lattice parameter, magnetic moment, band gap, and density of states have been simulated. U = 4.0 eV is recommended for AFM UO 2 , as this gives data close to experiments for all considered properties. For NM and AFM PuO 2 , we recommend U = 4.5 and 4.0 eV, respectively, though note that extremely large U values (c. 10 eV) are required to yield the most recently reported PuO 2 band gap. Exploration of the energies, lattice parameters, and band gaps of AFM UO 2 and NM PuO 2 calculated with PBESol + U (4.0 and 4.5 eV, respectively) + OMC reveals that several excited states have similar properties to the ground state and hence it is very hard to distinguish them from the ground state in the absence of OMC.